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Abstract. The aim of this work is to develop general optimization methods 
for finite difference schemes used to approximate linear differential equations. 
The specific case of the transport equation is exposed. In particular, the mini- 
mization of the numerical error is taken into account. The theoretical study of 
a related linear algebraic problem gives general results which can lead to the 
determination of the optimal scheme. 



1. Introduction: Scheme classes. Finite difference schemes used to approxi- 
mate linear differential equations induce numerical errors, that are generally diffi- 
cult to predict. The usual process consists in testing various schemes for more and 
' more refined time and space steps. 

. We here propose a completely different approach, which consists in determining the 

minimum norm error of a given finite difference scheme. This process has the ad- 
vantage of avoiding scheme convergence tests. Moreover, it can explain error jumps 

' that often occur in such approximations. 

O 

pH , Consider the transport equation: 

"tr« du du 

with the initial condition u(x,t = 0) = Uq(x). 
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A finite difference scheme for this equation can be written under the form: 

(2) 
where: 

ui m = u(lh,m,T) (3) 
I 6 {i — 1, i, i + 1}, m £ {n — 1, n, n + 1}, j ~ 0, n x , n ~ 0, n t , h, r 
denoting respectively the mesh size and time step. 
The Courant-Friedrichs-Lewy number (cfl) is defined as a = cr/h . 



A numerical scheme is specified by selecting appropriate values of the coefficients 
a, /3, 7, 5, e, Ci V> & an d $ in equation J5J. Values corresponding to numerical 
schemes retained for the present works are given in Tabled 
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Table 1. Numerical scheme coefficient. 
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The number of time steps will be denoted nt, the number of space steps, n x . In 
general, nt n x . 



The paper is organized as follows. The equivalent matrix equation is exposed in 
section Scheme optimization is presented in section [3] 



2. The Sylvester equation. 

2.1. Matricial form of the finite differences problem. Let us introduce the 
rectangular matrix defined by: 



U = [ Ui n ] !< 



i<7i x — l. l<n<nt 



(4) 



The problem JSJ) can be written under the following matricial form: 

MxU + U M 2 + £{U) = M 



(•5) 



where Mi, M2 and Mo are square matrices respectively n x — 1 by n x — 1, n t by nt, 
given by: 
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(7) X 



and where £ is a linear matricial operator which can be written as: 

£ = d + £2 + £3 + A 



(8) 
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where C\, £2, £3 and £4 are given by: 
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(9) 



(10) 



The second member matrix Mq bears the initial conditions, given for the specific 
value n = 0, which correspond to the initialization process when computing loops, 
and the boundary conditions, given for the specific values i = 0, i = n x . 

Denote by u exact the exact solution of 
The U exac t corresponding matrix will be: 



U t , 



where: 



with Xi = i h, t n = n r. 



\U exact i ] 0<i<n x — 1, 0<n<rif 



Uexacti — U eX act{^i-)^"n) 



(ii) 



(12) 



U is then solution of: 

Mi U + UM 2 + C(U) = M Q 
We will call error matrix the matrix defined by: 

E = U — U exac t 

Let us consider the matrix F defined by: 

F = Mi Uexact + Uexact M 2 + C(U eX act) ~ Mq 

The error matrix E satisfies then: 

Mi E + E M 2 + C{E) = F 



(13) 
(14) 

(15) 
(16) 



2.2. The matrix equation. 
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2.2.1. Theoretical formulation. Minimizing the error due to the approximation in- 
duced by the numerical scheme is equivalent to minimizing the norm of the matrices 
E satisfying ljiU)l . 

Since the linear matricial operator C appears only in the Crank-Nicholson scheme, 
we will restrain our study to the case C = 0. The generalization to the case C ^ 
can be easily deduced. 

The problem is then the determination of the minimum norm solution of: 



M 1 E + EM 2 = F 
which is a specific form of the Sylvester equation: 



(17) 



AX + XB = C (18) 

where A and B are respectively rabym and n by n matrices, C and X, m by n 
matrices. 

The solving of the Sylvester equation is generally based on Schur decomposition: 
for a given square n by n matrix A, n being an even number of the form n = 2p, 
there exists a unitary matrix U and a upper triangular block matrix T such that: 

A = U*TU (19) 

where U* denotes the (complex) conjugate matrix of the transposed matrix T U. 
The diagonal blocks of the matrix T correspond to the complex eigenvalues A; of 
A: 
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(20) 



T p ) 



where the block matrices Tj, i = 1, p are given by: 

Ke [Xi] Jm [Xi] 
- 1m [Xi] IZe [Xi] 

IZe being the real part of a complex number, and Im the imaginary one. 

Due to this decomposition, the Sylvester equation require, to be solved, that the 



(21) 



dimensions of the matrices be even numbers. We will therefore, in the following, 
restrain our study to n x and rit being even numbers. So far, it is interesting to 
note that the Schur decomposition being more stable for higher order matrices, it 
perfectly fits finite differences problems. 

Complete parametric solutions of the generalized Sylvester equation (??) is given 

As for the determination of the solution Sylvester equation, it is a major topic in 
control theory, and has been the subject of numerous works (see 01 [HI: El- 

na, mi, ca). 

In pp, the method is based on the reduction of the he observable pair (A,C) to 
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an observer-Hessenberg pair (H,D), H being a block upper Hessenberg matrix. 
The reduction to the observer-Hessenberg form (H, D) is achieved by means of the 
staircase algorithm (see @], ...). 

In |5] , in the specific case of B being a companion form matrix, the authors propose 
a very neat general complete parametric solution, which is expressed in terms of 
the controllability of the matrix pair (A,B), a symmetric matrix operator, and a 
parametric matrix in the Hankel form. 

We recall that a companion form, or Frobenius matrix is one of the following kind: 



/ o 
1 
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(22) 



These results can be generalized through matrix block decomposition to a block 
companion form matrix, which happens to be the case of our matrix M% in the 
specific case of n x and n t being even numbers: 

'' \ 

() M? 2 ... 



M 2 







V o 










(23) 



the M^ v , 1 < p < k being companion form matrices. 



Another method is presented in |14j. where the determination of the minimum- 
norm solution of a Sylvester equation is specifically developed. 
The accuracy and computational stability of the solutions is examined in [5] . 

2.2.2. Existence condition of the solution. Equation l|18ll has a unique solution if 
and only if A and B have no common eigenvalues. 

In our case, since M% is a upper triangular matrix whose diagonal coefficients are 
all equal to a, its eigenvalues are also all equal to a. As for the matrix Mi, one 
can easily check that a does not belong to its spectra. Hence, (|16f) has a unique 
solution, which accounts for the consistency of the given problem. 
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3. Scheme optimization. Advect a sinusoidal signal 

u = Cos [ ^ (a? -ct)] (24) 
A 

through the Lax scheme, where: 

A = n\ dx (25) 
n\ denotes the number of cells per wavelength. 

Let n\ remain unknown. 

Equation (|17|) can thus be normalized as: 

W[E + EM^ = F (26) 

where 

( m~i = 

{ M2 = 4^Af 2 (27) 

p _ hcfl p 

We deliberately choose a small value for the number of steps: n t = n x = 20, starting 
from the point that if the error is minimized for a small value of this number, it will 
be the same as this number increases. 

Figure displays the L2 norm of the isovalues of the error as a function of n\ (max- 
imums are in white, minimums in black; larger values are shown lighter). 

h 
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Figure 1 . Isovalues of the L2 norm of the error as a function of 
the number of cells per wavelength n\ 
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The square value of the L2 norm of the error, for two significative values of the 
number of cells per wavelength n p , is displayed in Figure [3 




Case 1 
Case 2 



20 40 60 80 100 — 

n 



Figure 2 . square value of the L 2 norm of the error as a function 
of the number of cells per wavelength n\. Case 1: A = 9. Case 2: 
A = 9.8. 

The above results ensure the faster convergence of the error. 

4. Conclusion. Thanks to the above results, we presently propose to optimize 
finite difference problems through minimization of the symbolic expression of the 
error as a function of the scheme parameters. 
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